### Jonathan D. Klingler, Gary E. Hollibaugh, Jr., and Adam J. Ramey
### "Don't Know What You Got: A Bayesian Hierarchical Model of Neuroticism and Ideological Uncertainty
### Political Science Research and Methods
###
###
### Main Analyses (All Figures except Figure 5; All Tables except B-1, B-2, and B-3)
### Description: This file creates (nearly) all tables and figures in the main paper.
###
### Note 1: Run the Preprocessing_and_Summary.R file beforehand to generate the CCES2014.Rda file.
### Note 2: Run the SingleCoreX.R files beforehand to generate the hierarchical model estimates.
### Note 3: Also run AMBoot.R beforehand to generate the bootstrapped Aldrich-McKelvey estimates.
### Note 4: Make sure to change the directory to the one containing this file.


rm(list=ls())
set.seed(1) # replicability


library(grDevices)
library(reshape2)
library(ggplot2)
library(RColorBrewer)
library(rjags)
library(plyr)
library(coda)
library(ggmcmc)
library(AER)
library(basicspace)
library(texreg)
library(grid)
library(Cairo)
library(xtable)

# load CCES data --- make sure in the correct directory first!
load("CCES2014.Rda")
cces_all <- cces
cces <- cces[!is.na(cces$self_emoti) & !is.na(cces$self_consc) & !is.na(cces$self_agree) 
             & !is.na(cces$self_extra) & !is.na(cces$self_openn) & !is.na(cces$CC421a) &
               !is.na(cces$newsint),]


# candidate names
candNames <- c("Obama","Cruz","Clinton","Paul","Bush","Democrats",
               "Republicans","Tea Party","Supreme Court")


# not sure and skipped are observationally (and arguably substantively)
# equivalent to each other
cces$totalNA <- as.numeric(cces$CC334C == "Not sure") + as.numeric(cces$CC334C == "Skipped") +  # obama
                as.numeric(cces$CC334D == "Not sure") + as.numeric(cces$CC334D == "Skipped") +  # clinton
                as.numeric(cces$CC334E == "Not sure") + as.numeric(cces$CC334E == "Skipped") +  # cruz
                as.numeric(cces$CC334F == "Not sure") + as.numeric(cces$CC334F == "Skipped") +  # rand paul
                as.numeric(cces$CC334G == "Not sure") + as.numeric(cces$CC334G == "Skipped") +  # jeb bush
                as.numeric(cces$CC334K == "Not sure") + as.numeric(cces$CC334K == "Skipped") +  # democratic party
                as.numeric(cces$CC334L == "Not sure") + as.numeric(cces$CC334L == "Skipped") +  # republican party
                as.numeric(cces$CC334M == "Not sure") + as.numeric(cces$CC334M == "Skipped") +  # tea party
                as.numeric(cces$CC334W == "Not sure") + as.numeric(cces$CC334W == "Skipped")    # scotus

# need to account for those who weren't asked about Cruz/Paul
cces$totalAsked <- 9 - (as.numeric(cces$CC334C == "Not Asked") + # obama
                        as.numeric(cces$CC334D == "Not Asked") + # clinton
                        as.numeric(cces$CC334E == "Not Asked") + # cruz
                        as.numeric(cces$CC334F == "Not Asked") + # rand paul
                        as.numeric(cces$CC334G == "Not Asked") + # jeb bush
                        as.numeric(cces$CC334K == "Not Asked") + # democratic party
                        as.numeric(cces$CC334L == "Not Asked") + # republican party
                        as.numeric(cces$CC334M == "Not Asked") + # tea party
                        as.numeric(cces$CC334W == "Not Asked")) # scotus


# rescaling the B5 to be between 0 and 1
cces$self_openn <- (cces$self_openn - 1)/6
cces$self_consc <- (cces$self_consc - 1)/6
cces$self_extra <- (cces$self_extra - 1)/6
cces$self_agree <- (cces$self_agree - 1)/6
cces$self_emoti <- (cces$self_emoti - 1)/6


# making the 150k or more level consistent
cces$faminc_recode <- cces$faminc
cces$faminc_recode[which(cces$faminc_recode == "$150,000 - $199,999")] <- levels(cces$faminc_recode)[17] 
cces$faminc_recode[which(cces$faminc_recode == "$200,000 - $249,999")] <- levels(cces$faminc_recode)[17] 
cces$faminc_recode[which(cces$faminc_recode == "$250,000 - $349,999")] <- levels(cces$faminc_recode)[17] 
cces$faminc_recode[which(cces$faminc_recode == "$350,000 - $499,999")] <- levels(cces$faminc_recode)[17] 
cces$faminc_recode[which(cces$faminc_recode == "$500,000 or more")] <- levels(cces$faminc_recode)[17] 
cces$faminc_recode[which(cces$faminc_recode == "$250,000 or more ")] <- levels(cces$faminc_recode)[17] 
cces$faminc_recode <- factor(cces$faminc_recode)


# total number of NAs                  
total.na.mod.1 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
                                                           self_agree + self_extra + self_openn, 
                                                           data = cces, family=binomial(link="probit"))
                                                           
total.na.mod.2 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
                                                           self_agree + self_extra + self_openn + 
                                                           I(gender == "Female") + I(2014 - birthyr) + 
                                                           I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                                           I(race == "Hispanic") + 
                                                           I(race == "Asian" | race == "Native American" | 
                                                             race == "Mixed" | race == "Other" | race == "Middle Eastern"), 
                                                           data = cces, family=binomial(link="probit"))
                                                           
total.na.mod.3 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
                                                           self_agree + self_extra + self_openn + 
                                                           I(gender == "Female") + I(2014 - birthyr) + 
                                                           I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                                           I(race == "Hispanic") + 
                                                           I(race == "Asian" | race == "Native American" | 
                                                             race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                                           as.numeric(educ), data = cces, family=binomial(link="probit"))
                                                           
total.na.mod.4 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
                                                           self_agree + self_extra + self_openn + 
                                                           I(gender == "Female") + I(2014 - birthyr) + 
                                                           I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                                           I(race == "Hispanic") + 
                                                           I(race == "Asian" | race == "Native American" | 
                                                             race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                                           as.numeric(educ) + I(newsint == "Most of the time") + 
                                                           I(newsint == "Don't know"), 
                                                           data = cces, family=binomial(link="probit"))
                                                          
total.na.mod.5 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
                                                           self_agree + self_extra + self_openn + 
                                                           I(gender == "Female") + I(2014 - birthyr) + 
                                                           I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                                           I(race == "Hispanic") + 
                                                           I(race == "Asian" | race == "Native American" | 
                                                             race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                                           as.numeric(educ) + I(newsint == "Most of the time") + 
                                                           I(newsint == "Don't know") + as.numeric(faminc_recode) + 
                                                           I(faminc_recode == "Prefer not to say"), 
                                                           data = cces, family=binomial(link="probit"))
                                                           
total.na.mod.6 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
                                                           self_agree + self_extra + self_openn + 
                                                           I(gender == "Female") + I(2014 - birthyr) + 
                                                           I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                                           I(race == "Hispanic") + 
                                                           I(race == "Asian" | race == "Native American" | 
                                                             race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                                           as.numeric(educ) + I(newsint == "Most of the time") + 
                                                           I(newsint == "Don't know") + as.numeric(faminc_recode) + 
                                                           I(faminc_recode == "Prefer not to say") + 
                                                           I(employ == "Full-time") + I(employ == "Part-time") + 
                                                           I(employ == "Unemployed") + I(employ == "Retired"), 
                                                           data = cces, family=binomial(link="probit"))

# coercing the models into a list for texreg
binom.mod.list <- list(total.na.mod.1, 
                       total.na.mod.2,
                       total.na.mod.3,
                       total.na.mod.4,
                       total.na.mod.5,
                       total.na.mod.6)
    
## table 1
# texreg for latex (commented out since it's already in the paper)                        
# texreg(binom.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
                          # custom.coef.names = c("Constant",
                                                # "Neuroticism",
                                                # "Conscientiousness",
                                                # "Agreeableness",
                                                # "Extraversion",
                                                # "Openness",
                                                # "Female",
                                                # "Age",
                                                # "Age^2/100",
                                                # "Black",
                                                # "Hispanic",
                                                # "Other Race",
                                                # "Education (1 = No HS; 6 = Postgrad)",
                                                # "High News Interest",
                                                # "Unknown News Interest",
                                                # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
                                                # "Income Refused",
                                                # "Employed Full-Time",
                                                # "Employed Part-Time",
                                                # "Unemployed",
                                                # "Retired"),
                          # reorder.coef = c(2:21,1),
                          # include.aic = FALSE, 
                          # include.deviance = FALSE)

# screenreg for review
# screenreg(binom.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
                          # custom.coef.names = c("Constant",
                                                # "Neuroticism",
                                                # "Conscientiousness",
                                                # "Agreeableness",
                                                # "Extraversion",
                                                # "Openness",
                                                # "Female",
                                                # "Age",
                                                # "Age^2/100",
                                                # "Black",
                                                # "Hispanic",
                                                # "Other Race",
                                                # "Education (1 = No HS; 6 = Postgrad)",
                                                # "High News Interest",
                                                # "Unknown News Interest",
                                                # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
                                                # "Income Refused",
                                                # "Employed Full-Time",
                                                # "Employed Part-Time",
                                                # "Unemployed",
                                                # "Retired"),
                          # reorder.coef = c(2:21,1),
                          # include.aic = FALSE, 
                          # include.deviance = FALSE)

# htmlreg for file
htmlreg(binom.mod.list, file = "table1.html", digits = 3, stars = c(0.01, 0.05, 0.1),
                          custom.coef.names = c("Constant",
                                                "Neuroticism",
                                                "Conscientiousness",
                                                "Agreeableness",
                                                "Extraversion",
                                                "Openness",
                                                "Female",
                                                "Age",
                                                "Age^2/100",
                                                "Black",
                                                "Hispanic",
                                                "Other Race",
                                                "Education (1 = No HS; 6 = Postgrad)",
                                                "High News Interest",
                                                "Unknown News Interest",
                                                "Income (1 = <10k; 12 = >150k; 13 = Refused)",
                                                "Income Refused",
                                                "Employed Full-Time",
                                                "Employed Part-Time",
                                                "Unemployed",
                                                "Retired"),
                          reorder.coef = c(2:21,1),
                          include.aic = FALSE, 
                          include.deviance = FALSE)


# calculating the model with the smallest BIC for plotting purposes
binom.min.bic <- which(c(BIC(total.na.mod.1), BIC(total.na.mod.2), 
                         BIC(total.na.mod.3), BIC(total.na.mod.4), 
                         BIC(total.na.mod.5), BIC(total.na.mod.6)) == 
                       min(c(BIC(total.na.mod.1), BIC(total.na.mod.2), 
                             BIC(total.na.mod.3), BIC(total.na.mod.4), 
                             BIC(total.na.mod.5), BIC(total.na.mod.6))))

# creating a mode function for the predicted probabilities
Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

## predicted probabilities for plotting (plot to be called later in combined plot)
# "some of the time" is the median value, and a plurality of people do NOT choose "most of the time"
# setting up the covariate values for prediction
total.na.preds <- predict(binom.mod.list[[binom.min.bic]], 
                          newdata = data.frame(self_emoti = seq(mean(cces$self_emoti) - 2*sd(cces$self_emoti), 
                                                                mean(cces$self_emoti) + 2*sd(cces$self_emoti), 
                                                                length.out = 100),
                                               self_consc = mean(cces$self_consc),
                                               self_agree = mean(cces$self_agree),
                                               self_extra = mean(cces$self_extra),
                                               self_openn = mean(cces$self_openn),
                                               newsint = levels(cces$newsint)[median(as.numeric(cces$newsint))],
                                               educ = Mode(cces$educ),
                                               employ = Mode(cces$employ),
                                               faminc_recode = Mode(cces$faminc_recode),
                                               gender = Mode(cces$gender),
                                               birthyr = median(cces$birthyr),
                                               race = Mode(cces$race)), se = TRUE)

# predicted probabilities with 90% CIs
# reversed scale because moving from emotional stability to neuroticism
total.na.predvals <- data.frame(neuroticism = seq(2, -2, length.out = 100),
                                low = pnorm(total.na.preds$fit - qnorm(0.95)*total.na.preds$se),
                                mid = pnorm(total.na.preds$fit),
                                high = pnorm(total.na.preds$fit + qnorm(0.95)*total.na.preds$se))


# standard aldrich-mckelvey
# need to remove those where a question wasn't asked (coded as 10)
# otherwise will make paul and cruz look too far right
newcces <- cces[!is.na(cces$self_emoti) & !is.na(cces$self_consc) & !is.na(cces$self_agree) 
                & !is.na(cces$self_extra) & !is.na(cces$self_openn) & !is.na(cces$CC421a) 
                & !is.na(cces$newsint) 
                & !(cces$CC334C == "Not Asked") 
                & !(cces$CC334D == "Not Asked") 
                & !(cces$CC334E == "Not Asked") 
                & !(cces$CC334F == "Not Asked") 
                & !(cces$CC334G == "Not Asked") 
                & !(cces$CC334K == "Not Asked") 
                & !(cces$CC334L == "Not Asked") 
                & !(cces$CC334M == "Not Asked") 
                & !(cces$CC334W == "Not Asked"),]
                
# pulling out the stimuli
aldmckspace <- cbind(as.numeric(newcces$CC334A), # self
                     as.numeric(newcces$CC334C), # obama
                     as.numeric(newcces$CC334D), # clinton
                     as.numeric(newcces$CC334E), # cruz
                     as.numeric(newcces$CC334F), # paul
                     as.numeric(newcces$CC334G), # bush
                     as.numeric(newcces$CC334K), # demparty
                     as.numeric(newcces$CC334L), # repparty
                     as.numeric(newcces$CC334M), # teaparty
                     as.numeric(newcces$CC334W)) # scotus

# naming the columns/stimuli
colnames(aldmckspace) <- c("Self", "Obama", "Clinton", "Cruz", "Paul", "Bush", "Democrats", "Republicans", "Tea Party", "Supreme Court")

# estimating the aldrich-mckelvey model
aldmck.out <- aldmck(aldmckspace, respondent = 1, polarity = 2, missing = c(8, 9))

# extracting the estimated "political information" from the A-M results
newcces$aldmck_polinfo <- aldmck.out$respondents$polinfo # aldrich-mckelvey political information (-1,1 correlation)


# tobit models with dependent variable as correlation (political information)
tobit.mod.1 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
                                      self_agree + self_extra + self_openn, 
                                      data = newcces, left = -1, right = 1)
                                      
tobit.mod.2 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
                                      self_agree + self_extra + self_openn + 
                                      I(gender == "Female") + I(2014 - birthyr) + 
                                      I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                      I(race == "Hispanic") + 
                                      I(race == "Asian" | race == "Native American" | 
                                        race == "Mixed" | race == "Other" | race == "Middle Eastern"), 
                                      data = newcces, left = -1, right = 1)
                                      
tobit.mod.3 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
                                      self_agree + self_extra + self_openn + 
                                      I(gender == "Female") + I(2014 - birthyr) + 
                                      I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                      I(race == "Hispanic") + 
                                      I(race == "Asian" | race == "Native American" | 
                                        race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                      as.numeric(educ), data = newcces, left = -1, right = 1)
                                      
tobit.mod.4 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
                                      self_agree + self_extra + self_openn + 
                                      I(gender == "Female") + I(2014 - birthyr) + 
                                      I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                      I(race == "Hispanic") + 
                                      I(race == "Asian" | race == "Native American" | 
                                        race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                      as.numeric(educ) + I(newsint == "Most of the time") + 
                                      I(newsint == "Don't know"), data = newcces, left = -1, right = 1)
                                      
tobit.mod.5 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
                                      self_agree + self_extra + self_openn + 
                                      I(gender == "Female") + I(2014 - birthyr) + 
                                      I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                      I(race == "Hispanic") + 
                                      I(race == "Asian" | race == "Native American" | 
                                        race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                      as.numeric(educ) + I(newsint == "Most of the time") + 
                                      I(newsint == "Don't know") + as.numeric(faminc_recode) + 
                                      I(faminc_recode == "Prefer not to say"), data = newcces, left = -1, right = 1)
                                      
tobit.mod.6 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
                                      self_agree + self_extra + self_openn + 
                                      I(gender == "Female") + I(2014 - birthyr) + 
                                      I((2014 - birthyr)^2/100) + I(race == "Black") + 
                                      I(race == "Hispanic") + 
                                      I(race == "Asian" | race == "Native American" | 
                                        race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
                                      as.numeric(educ) + I(newsint == "Most of the time") + 
                                      I(newsint == "Don't know") + as.numeric(faminc_recode) + 
                                      I(faminc_recode == "Prefer not to say") + 
                                      I(employ == "Full-time") + I(employ == "Part-time") + 
                                      I(employ == "Unemployed") + I(employ == "Retired"), 
                                      data = newcces, left = -1, right = 1)

# coercing the tobit models into a list for texreg
tobit.mod.list <- list(tobit.mod.1, 
                       tobit.mod.2,
                       tobit.mod.3,
                       tobit.mod.4,
                       tobit.mod.5,
                       tobit.mod.6)

## table 2
# texreg for latex (commented out since it's already in the paper)
# texreg(tobit.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
                          # custom.coef.names = c("Constant",
                                                # "Neuroticism",
                                                # "Conscientiousness",
                                                # "Agreeableness",
                                                # "Extraversion",
                                                # "Openness",
                                                # "Ln(scale)",
                                                # "Female",
                                                # "Age",
                                                # "Age^2/100",
                                                # "Black",
                                                # "Hispanic",
                                                # "Other Race",
                                                # "Education (1 = No HS; 6 = Postgrad)",
                                                # "High News Interest",
                                                # "Unknown News Interest",
                                                # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
                                                # "Income Refused",
                                                # "Employed Full-Time",
                                                # "Employed Part-Time",
                                                # "Unemployed",
                                                # "Retired"),
                          # reorder.coef = c(2:6, 8:22,1,7),
                          # include.aic = FALSE, 
                          # include.deviance = FALSE)

# screenreg for review
# screenreg(tobit.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
                          # custom.coef.names = c("Constant",
                                                # "Neuroticism",
                                                # "Conscientiousness",
                                                # "Agreeableness",
                                                # "Extraversion",
                                                # "Openness",
                                                # "Ln(scale)",
                                                # "Female",
                                                # "Age",
                                                # "Age^2/100",
                                                # "Black",
                                                # "Hispanic",
                                                # "Other Race",
                                                # "Education (1 = No HS; 6 = Postgrad)",
                                                # "High News Interest",
                                                # "Unknown News Interest",
                                                # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
                                                # "Income Refused",
                                                # "Employed Full-Time",
                                                # "Employed Part-Time",
                                                # "Unemployed",
                                                # "Retired"),
                          # reorder.coef = c(2:6, 8:22,1,7),
                          # include.aic = FALSE, 
                          # include.deviance = FALSE,
                          # include.wald = FALSE)


# html table
htmlreg(tobit.mod.list, file = "table2.html", digits = 3, stars = c(0.01, 0.05, 0.1),
                          custom.coef.names = c("Constant",
                                                "Neuroticism",
                                                "Conscientiousness",
                                                "Agreeableness",
                                                "Extraversion",
                                                "Openness",
                                                "Ln(scale)",
                                                "Female",
                                                "Age",
                                                "Age^2/100",
                                                "Black",
                                                "Hispanic",
                                                "Other Race",
                                                "Education (1 = No HS; 6 = Postgrad)",
                                                "High News Interest",
                                                "Unknown News Interest",
                                                "Income (1 = <10k; 12 = >150k; 13 = Refused)",
                                                "Income Refused",
                                                "Employed Full-Time",
                                                "Employed Part-Time",
                                                "Unemployed",
                                                "Retired"),
                          reorder.coef = c(2:6, 8:22,1,7),
                          include.aic = FALSE, 
                          include.deviance = FALSE,
                          include.wald = FALSE)


# finding the tobit model with the smallest BIC
tobit.min.bic <- which(c(BIC(tobit.mod.1), BIC(tobit.mod.2), 
                         BIC(tobit.mod.3), BIC(tobit.mod.4), 
                         BIC(tobit.mod.5), BIC(tobit.mod.6)) == 
                       min(c(BIC(tobit.mod.1), BIC(tobit.mod.2), 
                             BIC(tobit.mod.3), BIC(tobit.mod.4), 
                             BIC(tobit.mod.5), BIC(tobit.mod.6))))

# setting up the covariate values for prediction
tobit.preds <- predict(tobit.mod.list[[tobit.min.bic]], 
                       newdata = data.frame(self_emoti = seq(mean(newcces$self_emoti) - 2*sd(newcces$self_emoti), 
                                                             mean(newcces$self_emoti) + 2*sd(newcces$self_emoti), 
                                                             length.out = 100),
                                            self_consc = mean(newcces$self_consc),
                                            self_agree = mean(newcces$self_agree),
                                            self_extra = mean(newcces$self_extra),
                                            self_openn = mean(newcces$self_openn),
                                            newsint = levels(newcces$newsint)[median(as.numeric(newcces$newsint))],
                                            educ = Mode(newcces$educ),
                                            employ = Mode(newcces$employ),
                                            faminc_recode = Mode(newcces$faminc_recode),
                                            gender = Mode(newcces$gender),
                                            birthyr = median(newcces$birthyr),
                                            race = Mode(newcces$race)), se = TRUE, type = "linear")

# estimating the predicted correlations with 90% CIs
# reversed scale because moving from emotional stability to neuroticism
tobit.predvals <- data.frame(neuroticism = seq(2, -2, length.out = 100),
                                low = tobit.preds$fit - qnorm(0.95)*tobit.preds$se,
                                mid = tobit.preds$fit,
                                high = tobit.preds$fit + qnorm(0.95)*tobit.preds$se)


# combining the tobit plot with the binomial plot from before
binom.tobit.predvals <- rbind(data.frame(frame = total.na.predvals, plot = 1),
                              data.frame(frame = tobit.predvals, plot = 2))

# naming the columns of the data frame
colnames(binom.tobit.predvals) <- c("Neuroticism",
                                    "low",
                                    "mid",
                                    "high",
                                    "DV")

# setting up the appropriate result labels
binom.tobit.predvals$DV <- factor(binom.tobit.predvals$DV, 
                                  labels = c("Proportion of NA/DK Responses",
                                             "Correlation of Perceived Ideological Space with True Space"))

# figure 4
combined.plot <- ggplot(binom.tobit.predvals, aes(x = Neuroticism, y = mid, linetype = DV)) + 
                    geom_line(size = 0.75) + 
                    geom_ribbon(aes(ymin = low, ymax = high, fill = DV, linetype = NA), alpha = 0.6) + 
                    theme_bw() + 
                    xlab("Neuroticism\n(Number of Standard Deviations From Mean)") + 
                    ylab("Predicted Value of Dependent Variable\n") + 
                    scale_colour_discrete(name = c("Dependent Variable")) + 
                    scale_fill_grey(name = c("Dependent Variable"), start = 0.4, end = 0.8) +
                    scale_linetype_discrete(name = c("Dependent Variable")) + 
                    ylim(0,1) + 
                    theme(legend.position="bottom", legend.direction = "vertical")


cairo_ps(file="combined_plot.eps", height = 5, width = 7.5)
combined.plot
dev.off()




### load hierarchical model data
load("DKWYGCore1.Rda") 
load("DKWYGCore2.Rda") 
load("DKWYGCore3.Rda") 
load("DKWYGCore4.Rda") 
load("DKWYGCore5.Rda") 
load("DKWYGCore6.Rda")
load("DKWYGCore7.Rda") 
load("DKWYGCore8.Rda") 

thinning <- 4 # 12,500 per chain; 100,000 overall
a <- seq(1, dim(dkwyg.1[[1]])[1])
b <- a[seq(1, length(a), thinning)]


# coerce into a single object
DKWYG <- mcmc.list(
                   as.mcmc(dkwyg.1[[1]][b,]),
                   as.mcmc(dkwyg.2[[1]][b,]),
                   as.mcmc(dkwyg.3[[1]][b,]),
                   as.mcmc(dkwyg.4[[1]][b,]),
                   as.mcmc(dkwyg.5[[1]][b,]),
                   as.mcmc(dkwyg.6[[1]][b,]),
                   as.mcmc(dkwyg.7[[1]][b,]),
                   as.mcmc(dkwyg.8[[1]][b,])
                   )



# coerce the chains into a single chain for further analysis
BZ_data <- as.data.frame(runjags::combine.mcmc(DKWYG))

candNames <- c("Obama","Cruz","Clinton","Paul","Bush","Democrats",
               "Republicans","Tea Party","Supreme Court")
colnames(BZ_data)[86:94] <- candNames
melted_BZ <- melt(BZ_data)

# calculating posterior probability of parameter having same sign as posterior mean
postprob <- function(x){
   foo <- sum(sign(x) == sign(mean(x)))/length(x)
   return(foo)
   }

round(data.frame(apply(BZ_data, 2, postprob)), digits = 3)

   
# functions to create HPDs
middle95HPD <- function(x){
    out <- HPDinterval(as.mcmc(x), prob = 0.95)[,]
    names(out) <- c("ymin","ymax")
    return(out)
}

# alternative with middle 80%
middle80HPD <- function(x){
    out <- HPDinterval(as.mcmc(x), prob = 0.8)[,]
    names(out) <- c("ymin","ymax")
    return(out)
}

# alternative with middle 75%
middle75HPD <- function(x){
    out <- HPDinterval(as.mcmc(x), prob = 0.75)[,]
    names(out) <- c("ymin","ymax")
    return(out)
}


# alternative with middle 50%
middle50HPD <- function(x){
    out <- HPDinterval(as.mcmc(x), prob = 0.5)[,]
    names(out) <- c("ymin","ymax")
    return(out)
}



### distributional plots
# load bootstrapped AM replications
load("AMBoot.Rda") 

# rescale the replications so Obama = -1 and Cruz = 1
aldmck.boot.rescale <- 2*(aldmck.boot - aldmck.boot[,1])/(aldmck.boot[,3] - aldmck.boot[,1]) - 1

# use first 100,000 points since the thinned chains combined have 100,000 draws
aldmck.boot.rescale <- aldmck.boot.rescale[1:100000,]


# big plot for comparing distributions of estimates
melt.aldmck <- melt(aldmck.boot.rescale[,])
melt.aldmck <- melt.aldmck[,-1]
colnames(melt.aldmck)[1] <- "variable"
melt.aldmck$variable <- reorder(melt.aldmck$variable, melt.aldmck$value)
melted_BZ$variable <- reorder(melted_BZ$variable, melted_BZ$value)

# make sure order is correct!
aldmck.means <- aggregate(melt.aldmck[,2], list(melt.aldmck[,1]), mean)
aldmck.means <- aldmck.means[aldmck.means$Group.1 %in% candNames[-c(1,2)],]
aldmck.means <- aldmck.means[order(as.character(aldmck.means$Group.1)),]
BZ.means <- aggregate(melted_BZ[which(melted_BZ$variable %in% candNames),2], 
                      list(melted_BZ[which(melted_BZ$variable %in% candNames),1]), mean) 
BZ.means <- BZ.means[BZ.means$Group.1 %in% candNames[-c(1,2)],]
BZ.means <- BZ.means[order(as.character(BZ.means$Group.1)),]
r.squared <- summary(lm(aldmck.means[,2] ~ BZ.means[,2]))$r.squared

# combining the AM and BZ estimates into a data frame and labeling the methods accordingly
est_comparison_data <- rbind(data.frame(melt.aldmck, Method = "Aldrich-McKelvey"),
                             data.frame(melted_BZ[melted_BZ$variable %in% candNames,], Method = "Hierarchical Model"))

# figure 5
big_comparison_plot <- ggplot(est_comparison_data, aes(x=variable, y=value, fill = Method, shape = Method)) + 
                          stat_summary(data = est_comparison_data[est_comparison_data$variable %in% candNames[-c(1:2)],],
                                       fun.y=median,geom='point',position=position_dodge(width=0.725)) +
                          geom_violin(data = est_comparison_data[est_comparison_data$variable %in% candNames[-c(1:2)],], 
                                      aes(x= variable,y=value,fill=Method), alpha = 0.6, width = 0.75) + 
                          stat_summary(data = est_comparison_data[est_comparison_data$variable %in% candNames[-c(1:2)],],
                                       fun.y=median,geom='point',position=position_dodge(width=0.725)) +
                          stat_summary(data = est_comparison_data[est_comparison_data$variable %in% candNames[c(1:2)],],
                                       fun.y=median,geom='point') +
                          coord_flip() + 
                          theme_bw() +
                          theme(axis.title.y = element_blank()) +
                          scale_x_discrete(limits = rev(levels(factor(melt.aldmck[which(melt.aldmck$variable %in% candNames),]$variable)))) + 
                          scale_fill_manual(values = c("grey40", "white"), name = "Method") + 
                          ylab("\nEstimated Ideology") + 
                          annotate("text", x = 8, y = 1, label = paste("R^2 %~~%", sprintf("%.3f",round(r.squared, digits = 4))), parse = TRUE)


cairo_ps(file="bigcomparisonplot.eps", height = 5, width = 7.5)
big_comparison_plot
dev.off()




# examining the results of the hierarchical model
personality <- c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", 
              "Neuroticism")
personalityNews <- c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", 
              "Neuroticism","High News Interest")
varNames <- c("Constant", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", 
              "Neuroticism", "Female",
              "Age", "Age^2/100", "Black", "Hispanic", "Other Race",
              "Education", "High News Interest", "Unknown News Interest",
              "Income", "Income Refused",
              "Full Time Employment", "Part Time Employment", "Unemployed", "Retired")

colnames(BZ_data)[1:(length(varNames))] <- paste("beta_a[",varNames,"]",sep="")
colnames(BZ_data)[(length(varNames)+1):(2*length(varNames))] <- paste("beta_b[",varNames,"]",sep="")
colnames(BZ_data)[(2*length(varNames)+1):(3*length(varNames))] <- paste("beta_d[",varNames,"]",sep="")
colnames(BZ_data)[(3*length(varNames)+1):(4*length(varNames))] <- paste("beta_e[",varNames,"]",sep="")


# coerce into a separate object for testing
DKWYG.nofixed <- mcmc.list(
                           as.mcmc(dkwyg.1[[1]][b,-c(86,87)]),
                           as.mcmc(dkwyg.2[[1]][b,-c(86,87)]),
                           as.mcmc(dkwyg.3[[1]][b,-c(86,87)]),
                           as.mcmc(dkwyg.4[[1]][b,-c(86,87)]),
                           as.mcmc(dkwyg.5[[1]][b,-c(86,87)]),
                           as.mcmc(dkwyg.6[[1]][b,-c(86,87)]),
                           as.mcmc(dkwyg.7[[1]][b,-c(86,87)]),
                           as.mcmc(dkwyg.8[[1]][b,-c(86,87)])
                           )

gelman.diag(DKWYG.nofixed)

# print table a-1
table_a1_out <- gelman.diag(DKWYG.nofixed)
table_a1_psrf <- data.frame(table_a1_out$psrf[,2])
table_a1 <- data.frame(c(table_a1_psrf[[1]], table_a1_out$mpsrf))
rownames(table_a1) <- c(rownames(table_a1_psrf), "Multivariate PSRF")
rownames(table_a1)[1:(length(varNames))] <- paste("beta_a[",varNames,"]",sep="")
rownames(table_a1)[(length(varNames)+1):(2*length(varNames))] <- paste("beta_b[",varNames,"]",sep="")
rownames(table_a1)[(2*length(varNames)+1):(3*length(varNames))] <- paste("beta_d[",varNames,"]",sep="")
rownames(table_a1)[(3*length(varNames)+1):(4*length(varNames))] <- paste("beta_e[",varNames,"]",sep="")
rownames(table_a1) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a1)))))))
rownames(table_a1)[rownames(table_a1) == "c[5]"] <- "c[6]"
rownames(table_a1)[rownames(table_a1) == "c[4]"] <- "c[5]"
rownames(table_a1)[rownames(table_a1) == "c[3]"] <- "c[4]"
rownames(table_a1)[rownames(table_a1) == "c[2]"] <- "c[3]"
rownames(table_a1)[rownames(table_a1) == "c[1]"] <- "c[2]"
rownames(table_a1)[rownames(table_a1) == "gamma[3]"] <- "Clinton"
rownames(table_a1)[rownames(table_a1) == "gamma[4]"] <- "Paul"
rownames(table_a1)[rownames(table_a1) == "gamma[5]"] <- "Bush"
rownames(table_a1)[rownames(table_a1) == "gamma[6]"] <- "Democrats"
rownames(table_a1)[rownames(table_a1) == "gamma[7]"] <- "Republicans"
rownames(table_a1)[rownames(table_a1) == "gamma[8]"] <- "Tea Party"
rownames(table_a1)[rownames(table_a1) == "gamma[9]"] <- "Supreme Court"
colnames(table_a1) <- "PSRF Upper Bound"

print(xtable(table_a1), type = "html", file = "table-a1.html")



# hpd intervals
BZ_sum <- summary(as.mcmc(BZ_data))
hpd99 <- HPDinterval(as.mcmc(BZ_data), prob = 0.99)
hpd95 <- HPDinterval(as.mcmc(BZ_data), prob = 0.95)
hpd90 <- HPDinterval(as.mcmc(BZ_data), prob = 0.9)
hpd80 <- HPDinterval(as.mcmc(BZ_data), prob = 0.8)
hpd70 <- HPDinterval(as.mcmc(BZ_data), prob = 0.7)
hpd60 <- HPDinterval(as.mcmc(BZ_data), prob = 0.6)
hpd50 <- HPDinterval(as.mcmc(BZ_data), prob = 0.5)

hpd99.sig <- as.numeric(apply(hpd99,1, function(x){0 > max(x) | 0 < min(x)}))
hpd95.sig <- as.numeric(apply(hpd95,1, function(x){0 > max(x) | 0 < min(x)}))
hpd90.sig <- as.numeric(apply(hpd90,1, function(x){0 > max(x) | 0 < min(x)}))
hpd80.sig <- as.numeric(apply(hpd80,1, function(x){0 > max(x) | 0 < min(x)}))
hpd70.sig <- as.numeric(apply(hpd70,1, function(x){0 > max(x) | 0 < min(x)}))
hpd60.sig <- as.numeric(apply(hpd60,1, function(x){0 > max(x) | 0 < min(x)}))
hpd50.sig <- as.numeric(apply(hpd50,1, function(x){0 > max(x) | 0 < min(x)}))

smallestHPD <- apply(cbind(hpd99.sig, hpd95.sig, hpd90.sig, hpd80.sig, hpd70.sig, hpd60.sig, hpd50.sig),1,which.max)
neversig <- as.numeric(apply(cbind(hpd99.sig, hpd95.sig, hpd90.sig, hpd80.sig, hpd70.sig, hpd60.sig, hpd50.sig),1,sum) != 0)
smallestHPD[smallestHPD*neversig == 0] <- NA
smallestHPD <- as.numeric(gsub(".sig","",gsub("hpd","",colnames(cbind(hpd99.sig, hpd95.sig, hpd90.sig, hpd80.sig, hpd70.sig, hpd60.sig, hpd50.sig))[smallestHPD])))
smallestHPD.result <- data.frame(mean = BZ_sum$statistics[,'Mean'], se = BZ_sum$statistics[,'SD'], hpd80low = hpd80[,1], hpd80high = hpd80[,2], smallestHPD = smallestHPD)
rownames(smallestHPD.result) <- gsub("beta_s","sal_slo",gsub("beta_e","sal_int",gsub("beta_d","decisive",gsub("beta_b","op_slo",gsub("beta_a","op_int",rownames(smallestHPD.result))))))

# hpd intervals mentioned in text
round(HPDinterval(as.mcmc(BZ_data), prob = 0.69), digits = 3)
round(HPDinterval(as.mcmc(BZ_data), prob = 0.41), digits = 3)


# tables a-2 through a-8
table_a2 <- data.frame(hpd99)
rownames(table_a2) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a2)))))))
colnames(table_a2) <- c("Lower Bound", "Upper Bound")
rownames(table_a2)[rownames(table_a2) == "c[5]"] <- "c[6]"
rownames(table_a2)[rownames(table_a2) == "c[4]"] <- "c[5]"
rownames(table_a2)[rownames(table_a2) == "c[3]"] <- "c[4]"
rownames(table_a2)[rownames(table_a2) == "c[2]"] <- "c[3]"
rownames(table_a2)[rownames(table_a2) == "c[1]"] <- "c[2]"
print(xtable(table_a2, digits = 3), type = "html", file = "table-a2.html")

table_a3 <- data.frame(hpd95)
rownames(table_a3) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a3)))))))
colnames(table_a3) <- c("Lower Bound", "Upper Bound")
rownames(table_a3)[rownames(table_a3) == "c[5]"] <- "c[6]"
rownames(table_a3)[rownames(table_a3) == "c[4]"] <- "c[5]"
rownames(table_a3)[rownames(table_a3) == "c[3]"] <- "c[4]"
rownames(table_a3)[rownames(table_a3) == "c[2]"] <- "c[3]"
rownames(table_a3)[rownames(table_a3) == "c[1]"] <- "c[2]"
print(xtable(table_a3, digits = 3), type = "html", file = "table-a3.html")

table_a4 <- data.frame(hpd90)
rownames(table_a4) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a4)))))))
colnames(table_a4) <- c("Lower Bound", "Upper Bound")
rownames(table_a4)[rownames(table_a4) == "c[5]"] <- "c[6]"
rownames(table_a4)[rownames(table_a4) == "c[4]"] <- "c[5]"
rownames(table_a4)[rownames(table_a4) == "c[3]"] <- "c[4]"
rownames(table_a4)[rownames(table_a4) == "c[2]"] <- "c[3]"
rownames(table_a4)[rownames(table_a4) == "c[1]"] <- "c[2]"
print(xtable(table_a4, digits = 3), type = "html", file = "table-a4.html")

table_a5 <- data.frame(hpd80)
rownames(table_a5) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a5)))))))
colnames(table_a5) <- c("Lower Bound", "Upper Bound")
rownames(table_a5)[rownames(table_a5) == "c[5]"] <- "c[6]"
rownames(table_a5)[rownames(table_a5) == "c[4]"] <- "c[5]"
rownames(table_a5)[rownames(table_a5) == "c[3]"] <- "c[4]"
rownames(table_a5)[rownames(table_a5) == "c[2]"] <- "c[3]"
rownames(table_a5)[rownames(table_a5) == "c[1]"] <- "c[2]"
print(xtable(table_a5, digits = 3), type = "html", file = "table-a5.html")

table_a6 <- data.frame(hpd70)
rownames(table_a6) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a6)))))))
colnames(table_a6) <- c("Lower Bound", "Upper Bound")
rownames(table_a6)[rownames(table_a6) == "c[5]"] <- "c[6]"
rownames(table_a6)[rownames(table_a6) == "c[4]"] <- "c[5]"
rownames(table_a6)[rownames(table_a6) == "c[3]"] <- "c[4]"
rownames(table_a6)[rownames(table_a6) == "c[2]"] <- "c[3]"
rownames(table_a6)[rownames(table_a6) == "c[1]"] <- "c[2]"
print(xtable(table_a6, digits = 3), type = "html", file = "table-a6.html")

table_a7 <- data.frame(hpd60)
rownames(table_a7) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a7)))))))
colnames(table_a7) <- c("Lower Bound", "Upper Bound")
rownames(table_a7)[rownames(table_a7) == "c[5]"] <- "c[6]"
rownames(table_a7)[rownames(table_a7) == "c[4]"] <- "c[5]"
rownames(table_a7)[rownames(table_a7) == "c[3]"] <- "c[4]"
rownames(table_a7)[rownames(table_a7) == "c[2]"] <- "c[3]"
rownames(table_a7)[rownames(table_a7) == "c[1]"] <- "c[2]"
print(xtable(table_a7, digits = 3), type = "html", file = "table-a7.html")

table_a8 <- data.frame(hpd50)
rownames(table_a8) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a8)))))))
colnames(table_a8) <- c("Lower Bound", "Upper Bound")
rownames(table_a8)[rownames(table_a8) == "c[5]"] <- "c[6]"
rownames(table_a8)[rownames(table_a8) == "c[4]"] <- "c[5]"
rownames(table_a8)[rownames(table_a8) == "c[3]"] <- "c[4]"
rownames(table_a8)[rownames(table_a8) == "c[2]"] <- "c[3]"
rownames(table_a8)[rownames(table_a8) == "c[1]"] <- "c[2]"
print(xtable(table_a8, digits = 3), type = "html", file = "table-a8.html")




### table 3 information 
table3 <- as.data.frame(smallestHPD.result[!(rownames(smallestHPD.result) %in% candNames),-c(3,4)])
rownames(table3) <- gsub("tau", "c", gsub("sal_slo","Saliency Slope",gsub("sal_int","Saliency Intercept",gsub("decisive","Decisiveness",gsub("op_slo","Opinion Slope",gsub("op_int","Opinion Intercept",rownames(table3)))))))
rownames(table3)[rownames(table3) == "c[5]"] <- "c[6]"
rownames(table3)[rownames(table3) == "c[4]"] <- "c[5]"
rownames(table3)[rownames(table3) == "c[3]"] <- "c[4]"
rownames(table3)[rownames(table3) == "c[2]"] <- "c[3]"
rownames(table3)[rownames(table3) == "c[1]"] <- "c[2]"
table3 <- rbind(table3, NA) 
rownames(table3)[dim(table3)[1]] <- "Number of Observations"
table3[dim(table3)[1],1] <- dim(newcces)[1]
colnames(table3) <- c("Posterior Mean", "Posterior SD", "Best Nonzero HPD")
print(xtable(table3, digits = 3), type = "html", file = "table3.html")



melted_BZ_labeled <- melt(BZ_data)

#melted_data[grep("*[6]",melted_data$variable),]$value <- -melted_data[grep("*[6]",melted_data$variable),]$value

melted_BZ_labeled <- melted_BZ_labeled[(1:nrow(melted_BZ_labeled))[-grep("Constant", melted_BZ_labeled$variable)],]

# extract results dealing only with personality
melted_BZ_B5 <- melted_BZ_labeled[grep(paste(personality,collapse="|"), melted_BZ_labeled$variable),]
  


### one plot for gb (opinion slopes), ge (saliency intercepts), and gd (decisiveness)
bigplot.frame <- rbind(data.frame(value = melted_BZ_B5[grep("beta_b", melted_BZ_B5$variable),], class = 1),
                       data.frame(value = melted_BZ_B5[grep("beta_e", melted_BZ_B5$variable),], class = 2),
                       data.frame(value = melted_BZ_B5[grep("beta_d", melted_BZ_B5$variable),], class = 3))

colnames(bigplot.frame) <- c("variable", "value", "class")
bigplot.frame$class <- factor(bigplot.frame$class, labels = c("Opinion Slopes", "Saliency Intercepts", "Decisiveness"))
levels(bigplot.frame$variable) <- list(Openness = c("beta_b[Openness]", "beta_e[Openness]", "beta_d[Openness]"),
                                       Conscientiousness = c("beta_b[Conscientiousness]", "beta_e[Conscientiousness]", "beta_d[Conscientiousness]"),
                                       Extraversion = c("beta_b[Extraversion]", "beta_e[Extraversion]", "beta_d[Extraversion]"),
                                       Agreeableness = c("beta_b[Agreeableness]", "beta_e[Agreeableness]", "beta_d[Agreeableness]"),
                                       Neuroticism = c("beta_b[Neuroticism]", "beta_e[Neuroticism]", "beta_d[Neuroticism]"))

# figure 6
bigplot <- ggplot(bigplot.frame, aes(x=variable,y=value)) + 
    geom_violin(alpha=0.6) +
    theme_bw() +
    facet_grid(~ class) + 
    #  uncomment this line and comment out the one below for a plot with all variables
    #  scale_x_discrete(labels=varNames[-1]) +
    scale_x_discrete(labels=varNames[c(2:6,14)]) +
    theme(axis.title.y = element_blank(), legend.title=element_blank(), legend.position="none", panel.margin = unit(1, "lines")) +
    ylab("\nCoefficient Estimate") +
    geom_hline(yintercept = 0, linetype="dotted") +
    coord_flip() + 
    scale_fill_manual(values = "white") +
    stat_summary(fun.y=median,geom='point') + 
    stat_summary(fun.data=middle50HPD,geom='errorbar',width = 0.075) +
    stat_summary(fun.y=middle80HPD,geom='line',size = 0.5)

cairo_ps(file="bigplot.eps", height = 4, width = 8)
bigplot
dev.off()

# how much of the PD for the Neuroticism coefficient in the opinion slope equations lies to the left of zero?
sum(bigplot.frame[which(bigplot.frame$class == "Opinion Slopes" & bigplot.frame$variable == "Neuroticism"),2] < 0)/
length(bigplot.frame[which(bigplot.frame$class == "Opinion Slopes" & bigplot.frame$variable == "Neuroticism"),2] < 0)

# saliency intercept?
sum(bigplot.frame[which(bigplot.frame$class == "Saliency Intercepts" & bigplot.frame$variable == "Neuroticism"),2] < 0)/
length(bigplot.frame[which(bigplot.frame$class == "Saliency Intercepts" & bigplot.frame$variable == "Neuroticism"),2] < 0)

# and for Decisiveness?
sum(bigplot.frame[which(bigplot.frame$class == "Decisiveness" & bigplot.frame$variable == "Neuroticism"),2] < 0)/
length(bigplot.frame[which(bigplot.frame$class == "Decisiveness" & bigplot.frame$variable == "Neuroticism"),2] < 0)



### appendix plots (PNG used because EPS is way too slow)
plot.labs <- data.frame(Parameter = c(colnames(DKWYG.nofixed[[1]])), Label = c(rep(varNames,4),"Copartisan",candNames[-c(1:2)],"Cutpoint 2", "Cutpoint 3", "Cutpoint 4", "Cutpoint 5", "Cutpoint 6"))

# opinion intercept
CairoPNG(file = "opinion-intercept-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_a", burnin = F)) + 
              theme_bw() + 
              ylab('Value\n') + 
              xlab('\nIteration') + 
              facet_wrap(~Parameter, ncol = 3, scales = "free") + 
              geom_line(alpha = 0.05)
dev.off()


CairoPNG(file = "opinion-intercept-density.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_a", burnin = F)) + 
            theme_bw() + 
            ylab('Density\n') + 
            xlab('\nValue') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free")
dev.off()


CairoPNG(file = "opinion-intercept-running.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_a", burnin = F)) + 
            theme_bw() + 
            ylab('Value\n') + 
            xlab('\nIteration') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free") + 
            geom_line(alpha = 0.3)
dev.off()




# opinion slope   
CairoPNG(file = "opinion-slope-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)           
ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_b", burnin = F)) + 
              theme_bw() + 
              ylab('Value\n') + 
              xlab('\nIteration') + 
              facet_wrap(~Parameter, ncol = 3, scales = "free") + 
              geom_line(alpha = 0.05)
dev.off()

CairoPNG(file = "opinion-slope-density.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_b", burnin = F)) + 
            theme_bw() + 
            ylab('Density\n') + 
            xlab('\nValue') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free")
dev.off()

CairoPNG(file = "opinion-slope-running.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_b", burnin = F)) + 
            theme_bw() + 
            ylab('Value\n') + 
            xlab('\nIteration') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free") + 
            geom_line(alpha = 0.3)
dev.off()





# decisiveness slope    
CairoPNG(file = "decisiveness-slope-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                     
ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_d", burnin = F)) + 
              theme_bw() + 
              ylab('Value\n') + 
              xlab('\nIteration') + 
              facet_wrap(~Parameter, ncol = 3, scales = "free") + 
              geom_line(alpha = 0.05)
dev.off()              


CairoPNG(file = "decisiveness-slope-density.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_d", burnin = F)) + 
            theme_bw() + 
            ylab('Density\n') + 
            xlab('\nValue') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free")
dev.off()


CairoPNG(file = "decisiveness-slope-running.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_d", burnin = F)) + 
            theme_bw() + 
            ylab('Value\n') + 
            xlab('\nIteration') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free") + 
            geom_line(alpha = 0.3)
dev.off()




# saliency intercept    
CairoPNG(file = "saliency-intercept-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                               
ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_e", burnin = F)) + 
              theme_bw() + 
              ylab('Value\n') + 
              xlab('\nIteration') + 
              facet_wrap(~Parameter, ncol = 3, scales = "free") + 
              geom_line(alpha = 0.05)
dev.off()

CairoPNG(file = "saliency-intercept-density.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_e", burnin = F)) + 
            theme_bw() + 
            ylab('Density\n') + 
            xlab('\nValue') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free")
dev.off()


CairoPNG(file = "saliency-intercept-running.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_e", burnin = F)) + 
            theme_bw() + 
            ylab('Value\n') + 
            xlab('\nIteration') + 
            facet_wrap(~Parameter, ncol = 3, scales = "free") + 
            geom_line(alpha = 0.3)
dev.off()




# saliency slope     
CairoPNG(file = "saliency-slope-traceplot.png", height = 2.5, width = 7.5, units = "in", dpi = 300)                                        
ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_s", burnin = F)) + 
              theme_bw() + 
              ylab('Value\n') + 
              xlab('\nIteration') + 
              geom_line(alpha = 0.05) + 
            facet_wrap(~Parameter, scales = "free")
dev.off()              

CairoPNG(file = "saliency-slope-density.png", height = 2.5, width = 7.5, units = "in", dpi = 300)
ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_s", burnin = F)) + 
            theme_bw() + 
            ylab('Density\n') + 
            xlab('\nValue')  + 
            facet_wrap(~Parameter, scales = "free")
dev.off()              


CairoPNG(file = "saliency-slope-running.png", height = 2.5, width = 7.5, units = "in", dpi = 300)
ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_s", burnin = F)) + 
            theme_bw() + 
            ylab('Value\n') + 
            xlab('\nIteration') + 
            geom_line(alpha = 0.3) + 
            facet_wrap(~Parameter, scales = "free")
dev.off()



              
# candidates    
CairoPNG(file = "stimuli-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                                                  
ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "gamma", burnin = F)) + 
              theme_bw() + 
              ylab('Value\n') + 
              xlab('\nIteration') + 
              facet_wrap(~Parameter, ncol = 1, scales = "free") + 
              geom_line(alpha = 0.05)
dev.off()              
           
CairoPNG(file = "stimuli-density.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "gamma", burnin = F)) + 
            theme_bw() + 
            ylab('Density\n') + 
            xlab('\nValue') + 
            facet_wrap(~Parameter, ncol = 1, scales = "free")
dev.off()

CairoPNG(file = "stimuli-running.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "gamma", burnin = F)) + 
            theme_bw() + 
            ylab('Value\n') + 
            xlab('\nIteration') + 
            facet_wrap(~Parameter, ncol = 1, scales = "free") + 
            geom_line(alpha = 0.3)
dev.off()


              
# cutpoints  
CairoPNG(file = "cutpoints-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                                                              
ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "tau", burnin = F)) + 
              theme_bw() + 
              ylab('Value\n') + 
              xlab('\nIteration') + 
              facet_wrap(~Parameter, ncol = 1, scales = "free") + 
              geom_line(alpha = 0.3)
dev.off()              

CairoPNG(file = "cutpoints-density.png", height = 10, width = 10, units = "in", dpi = 300)              
ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "tau", burnin = F)) + 
            theme_bw() + 
            ylab('Density\n') + 
            xlab('\nValue') + 
            facet_wrap(~Parameter, ncol = 1, scales = "free")
dev.off()            

CairoPNG(file = "cutpoints-running.png", height = 10, width = 10, units = "in", dpi = 300)
ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "tau", burnin = F)) + 
            theme_bw() + 
            ylab('Value\n') + 
            xlab('\nIteration') + 
            facet_wrap(~Parameter, ncol = 1, scales = "free") +
            geom_line(alpha = 0.3)
dev.off()

  
